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SOLUTION OF ELASTOPLASTIC TORSION PROBLEM 
BY BOUNDARY INTEGRAL METHOD 
by Alexander Mendelson 
Lewis Research Center 

SUMMARY 

The boundary integral method was applied to the elastoplastic analysis of the torsion 
of prismatic bars. The boundary of the bar cross section was divided into n intervals, 
and the unknown warping function was assumed constant over each interval. Since the 
unknowns of the problem appear only on the boundary, the dimensions of the problem 
are reduced effectively by one. 

Numerical results were obtained for a bar with a square cross section, and com- 
parisons were made with results obtained by the finite difference method. Very good 
accuracy was obtained using relatively small sets of linear algebraic equations. The 
convergence rate for the elastoplastic problem was also much faster using the boundary 
integral method than using the finite difference method. This should be reflected in an" 
appreciable savings in computer time. The method can be readily programmed for a 
digital computer for both simply and multiply connected bodies. 


INTRODUCTION 

The problem of elastic torsion of prismatic bars is one of the classical problems of 
mechanics and as such has received extensive treatment. With the advent of the high 
speed digital computer the corresponding elastoplastic problem has also come to the 
forefront. 

The usual solutions to the elastoplastic torsion problem assume perfect plasticity 
and are very often limited to cases in which the complete cross section is plastic, be- 
cause the elastoplastic boundary is considered difficult to find. Recently, general solu- 
tions have been presented by the author (refs. 1 and 2) using finite differences and the 
method of elastic solutions and by Herakovich (refs. 3 and 4) using a minimum rate 
principle of plasticity in conjunction with nonlinear programming. In reference 4 



Herakovich compares his results with those of reference 2 and concludes that the mini- 
mum rate principle would require a great many more finite elements for its implemen- 
tation in order to approach the accuracy of the finite difference method. 

More recently a great deal of interest has been aroused by a new approach involving 
boundary integral techniques to solve elastic and elastoplastic problems. A general 
review of the boundary integral techniques is presented in reference 5. The first 
numerical application of these techniques to an elastoplastic problem involving a notched 
beam in bending is described in reference 6. A major advantage of these techniques is 
that a problem need only be solved for the unknowns on the boundary of the cross section 
of the body. This effectively reduces the dimensions of the problem by one. 

The development of the theory of the boundary integral method to the torsion prob- 
lem is described in general terms in reference 5. The purpose of the present report is 
to show specifically how the theory is applied to the elastoplastic torsion of bars of any 
geometry. Numerical computations are presented for bars having square cross sec- 
tions. Comparisons are made for the elastic case with the known analytical solution and 
for the elastoplastic case with the finite difference solution. 


METHOD OF ANALYSIS 

The elastoplastic torsion problem can be formulated in several ways as shown in 
reference 5. In particular, for example, it is shown in reference 5 that, by formulating 
the problem in terms of a stress function a closed form solution can be obtained for a 
circular bar with linear strain hardening. In general, however, it would seem that a 
formulation in terms of the warping function should be preferable, since the warping 
function (axial displacement) is physically more meaningful than the stress function and, 
more importantly, the distinction between simply connected and multiply connected re- 
gions disappears (ref. 7). The warping function formulation was therefore used herein. 


Basic Equations 

The basic equations necessary for the formulation to be used are as follows (see 
fig. 1 for the coordinate system): 

Equilibrium equation: 


dr. 


xz 


dr. 


dx 


I5 = o 
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where t and t, 7 „ are the shear stresses acting on the cross section. 
Stress -strain relations: 


T xz 2< ^ € xz ” e xz^ 


r\ 


T yz 2 ^- c yz " € yz^ 


J 


(2) 


where G is the shear modulus, e xz and e yz are the total (tensor) strains, and 

and €*L are the plastic shear strains. It should be noted that equation (2) implies using 
yz 

the deformation theory of plasticity. However, as shown in reference 8, both the total 
and incremental theories of plasticity furnish the same solution to the torsion problem 
provided either the cross section is circular or the material is perfectly plastic. It is 
reasonable to assume, therefore, that this will be approximately true for most practical 
problems. Indeed, it has been shown by Huth in reference 9 that for the case of a square 
cross section with strain hardening there is little difference between incremental and 
deformation theories. In what follows, therefore, use will be made of the deformation 
or total theories of plasticity, and the load will be assumed to be applied in one step. 

The use of incremental theories does not appreciably complicate the problem, and the 
necessary formulation is given in reference 2 for those desiring to use it. 

Saint- Venant relations: 




(3) 


where a is the angle of twist per unit length and w the warping function (the axial dis- 
placement). 

Substituting equation (3) into equation (2) and the resulting equation into equation (1) 
result in 


3?Z + i?2!, F ( X;y) 

dx 2 ay 2 


where 


F(x, y) = 2 



(4a) 


(4b) 
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The plastic strains appearing in the definition of the function F(x, y) are of course 
in turn nonlinear functions of warping function w. They can be determined from (ref. 1) 



€ P =_P € 
yz , yz 
e et J 


( 5 ) 


where 





( £ xz> 2 + 


(6a) 


and 


e p = fle et> < 6b > 

Equation (6b) represents the uniaxial stress- strain curve in terms of equivalent plastic 
strain against equivalent total strain; that is, 


e p = e et 


1% 
3 G 


(6c) 


where a g , the equivalent stress, represents the stress on the uniaxial stress-strain 
curve and e the plastic strain on that curve. Thus, for a given stress -strain curve, 

IT 

the relation between e p and represented by equation (6b) can be determined using 
equation (6c). For the case of linear strain hardening, the relation (6b) can be written 
as 


<P = 


e et - - (1 + M) — 
3 E_ 

i + 2(i + p)_a_ 
3 1 - m 


(6d) 



where a Q is the yield stress, ju Poisson's ratio, and m the linear strain hardening 
parameter (i.e. , the ratio of the slope of the linear strain hardening line to the elastic 
modulus). 


Boundary Conditions 

The boundary conditions for an unloaded lateral surface are given by 

* T xz + mT y Z = 0 (7) 

where l and m are the direction cosines of the external normal to the surface with 
respect to the x and y axes, respectively. Substituting equations (2) and (3) into equa- 
tion (7) and malting use of the relation 

7 3w m Sw 3w 
L — + m — = — 

3x 3y . 3n 


result in 


— = w' = a(ly - mx) + 2(Ze£ z + mej z ) (8) 

9n “ 

where w' is the outward normal derivative of the warping function w at the boundary of 
the surface. For a rectangular boundary parallel to one of the coordinate axes, the sec- 
ond term on the right side of equation (8) vanishes. This follows from the fact that if 

l - 0, then e)L = 0, and if m = 0, then = 0. 

y & xzi 


Boundary Integral Formulation 

The differential equation (4a) can be readily converted to a boundary integral equa- 
tion by using Green’s second theorem as shown in reference 5. The resultant equation 
is 


W<P) ^// 


F(Q) In r pQ dA + - 


- f w(q)-L 

77 1 0n q 


(In r p q ) dq 


;/ w ' 


(q) In r pq dq 


(9) 
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The notation in equation (9) is fully described in figure 2. The capital letters P and Q 
denote interior points while the lower case letters p and q denote boundary points. 
Thus, for example, r p Q represents the distance between the boundary point p and the 
interior point Q. The outward normal derivative at the boundary point q is 9/ 3n q . 

The same equation applies for a point P in the interior of the region R except that the 
right side of the equation is divided by 2 ; that is, 




w(P) = — / F(Q) In r pQ dA 


i/ 


w(q) — (In r p ) dq - — 
9n 2 


L [. 

’’Jr 


(q) In r pq dq 
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In equation (9), the normal derivative w'(q) is known from equation (8). If we 
assume that the function F defined by equation (4b) is also known, then equation (9) 
represents an integral equation for the unknown function w(q) on the boundary C of the 
region R. Once this integral equation is solved for w(q), the warping function w can 
be calculated at any interior point P by direct quadrature using equation (10). 

Actually, however, the function F(Q) is not known to start with, since it is a func- 
tion of the plastic strains which are determined as part of the solution. However, F(Q) 
can be determined by an iterative process. First, equation (9) is solved for an arbi- 
trary distribution (say zero) of F(Q). Then w(P) is calculated in the interior of the 
region by means of equation (10). The total strains are computed from equation (3), and 
the plastic strains are computed by means of equations (5) and (6). A better approxima- 
tion to the function F is then obtained from equation (4b). Equation (9) is now' re- 
solved, and the process is repeated until convergence is obtained. 


Numerical Procedure 

To solve equation (9) for the unknown function w(p), the straightforward procedure 
of replacing the integrals by summations can be used. The boundary is divided into n 
intervals with a nodal point taken at the center of each interval. The unknown function is 
assumed constant over each interval. Similarly, the region R is divided into a number 
of cells and the function F assumed constant over each cell. Equation (9) is then 
written for each nodal point as follows: 

£ <a ij - v )w j “ £ vi + r i 

j=i M 
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i = 1,2,. . .,n 


(ID 



where is the Kronecker delta (equal to 0 when i * j and equal to 1 when i = j). 

The coefficients a^, b^, and are given in the appendix. We thus have n equations 
for the n unknowns w-. This set of equations can readily be solved by any standard 
procedure. 

Once the w^ are known on the boundary, equation (10) can be used to calculate w 
at any interior point. However, as can be seen from equation (3), although w is not 
really needed to calculate the strains, its derivatives are. To determine the derivatives 
of w, equation (10) can be differentiated directly to give 




For 3w(P)/3y we interchange x and y. 

Again we replace the integrals by sums and write 


‘i Z 'VkV - »k m k B ijk - «k c iik> * Z F w D iju 

^Lk=l k, l 


where the coefficients A^ k , B^, C^, and are listed in the appendix, Z k and 

m k are the direction cosines of the normal at the k tl1 boundary interval, and £ is 
the sum for all the plastic cells in the region. k, l 

When the derivatives of w are used, the total strains are computed from equa- 
tion (3), the plastic strains from equation (5) and (6), and the function F from equa- 
tion (4b). The process is then repeated until convergence is obtained. The degree of 
convergence can be determined from a relation of the form 
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where M is the total number of points P i flowing plastically and k- 1 and k are two 
successive iterations. The convergence criterion K can be made as small as desired. 
For the calculations reported herein values of 10"^ and 10 ® were used for K. 

For a cross section symmetric about one or both of the coordinate axes, advantage 
can be taken of the symmetries (or antisymmetries) to reduce the number of unknowns. 
For example, for the square plate for which calculations are reported herein, symmetry 
exists not only about the x and y axes, but also about the diagonals. Therefore, it is 
necessary to consider only the boundary between A and B (fig. 3) since the values of w 
are the same along AB, CD, EF, and GH, and these are the negatives of the values 
along CB, ED, GF, and AH. Similarly, the function F, defined by equation (4b), is 
antisymmetric with respect to both the x and y axes. 


RESULTS AND DISCUSSION 

Calculations were performed by the previous technique for a bar of square cross 
section as shown in figure 3. The dimensionless angle of twist per unit length { 3 , de- 
fined as ora/e 0 (where e Q is the yield strain and a is 1/2 the side of the square), was 
increased in steps of one from j3 = 1 to jS = 6. Linear strain hardening was assumed 
with values of the strain hardening parameter taken as 0 (perfect plasticity), 0.05, 0. 1, 
and 0. 2. Poisson's ratio was assumed as 0. 3 in all calculations. 

For j3 = 1, the bar is elastic and a comparison was made between the analytical 
solution as given, for example, in reference 7 as well as with the finite difference solu- 
tion of reference 2. The results are shown in tables I to IH. Table I shows the warping 
function as computed on the boundary of the bar cross section. The comparison with the 
analytical solution of reference 7 shows very good agreement with just 4 unknowns to 
solve for in the boundary integral method. The analytical solution of reference 7 uses 
an infinite series of which 30 terms were summed. Summing 50 terms did not change 
the answers to 5 significant figures. The warping function was not computed in refer- 
ence 2. 

Table n shows the comparison for the maximum shear stress (at the center of the 
edge of the square) and the moment with the analytical solution of reference 7 and the 
finite difference solution of reference 2. Again it is seen that with just 4 unknowns in 
the boundary integral method very good results are obtained, as good as the results ob- 
tained for the finite difference method using 55 unknowns. 
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Table in presents the dimensionless shear stress distribution in the x-direction 
(r = r /2 GeJ throughout the cross section using 10 unknowns for the boundary integral 
method and 55 unknowns for the finite difference method. Again excellent agreement 
was obtained. Actually, the results with 4 unknowns using the boundary integral method 
are almost as good, but the results with 10 unknowns are presented to match the actual 
(x, y) values of the finite difference results without having to cross plot. 

The dimensionless angle of twist per unit length /3 was then increased in unit steps 
to a maximum value of /3 = 6 for each value of the strain hardening parameter m. The 
total boundary was divided into 80 intervals resulting in 10 equations for 10 unknowns. 
Several test calculations were made with fewer intervals, and the results indicated that 
using 48 intervals (6 unknowns) changed the moment and maximum stress by at most one 
in the third significant figure and changed the maximum plastic strain by about 3 percent. 
All the subsequent results are therefore shown for 80 intervals (10 unknowns), although 
from an engineering viewpoint 48 or even 32 intervals would be sufficient. 

The results of the calculations are summarized in table IV and figures 4 to 6. Fig- 

O 

ure 4 shows the dimensionless moment defined as M* = M/(2Ge 0 a ) for various values 
of /3 and m. Figure 5 shows the corresponding dimensionless maximum shear 
stresses defined as 


and figure 6 shows the spread of the plastic zones with an increase of the angle of 
twist 0. 

In all the calculations the convergence number K of equation (14) was taken either 
as 0. 0001 or 0. 00001. In many of the calculations both numbers were used in turn. The 
differences in the results were found to be insignificant. For example, the number of 
iterations for convergence for the case of maximum plastic flow, which occurred for 
(3 = 6 and m = 0, was 39 for K = 0. 0001 and 53 for K = 0. 00001, and the results were 
all the same to at least three significant figures. For the case (3 = 5 and m = 0, the 
number of iterations for K = 0. 0001 was 33. For the same case using finite differences, 
203 iterations were required. 


CONCLUDING REMARKS 

The boundary integral method was found to be very suitable for the elastoplastic 
analysis of the torsion of prismatic bars. Very good accuracy can be obtained by using 
relatively small sets of linear algebraic equations. 


9 



A comparison with the finite difference method indicates a great savings in the num- 
ber of unknowns that have to be determined and also a much faster convergence rate 
using the method of successive elastic solutions for both formulations. This should be 
reflected in appreciable savings in computer time. 

The boundary integral method can readily be programmed in a straightforward man- 
ner for a digital computer. The use of the warping function to formulate the problem 
permits applying the method with equal ease to both simply connected and multiply con- 
nected bodies . 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 30, 1974, 

505-01. 
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APPENDIX - BOUNDARY INTEGRAL COEFFICIENTS 


The division of the boundary into intervals with their corresponding nodal points is 
shown in figure 7. The x and y coordinates of a boundary nodal point are desig- 
nated as (x^., y bi ). The coordinates at the beginning and end of an interval (say interval 
j) are designated by (^, r?j) at the beginning of the interval and by (£. +1 , T?j + j) at the enc * 
of the interval. The interval lengths h. need not be equal. The coordinates of the cen- 
troid of an interior cell where plastic flow occurs are designated by (x k , x^). 

The coefficients in equation (11) are then given by 


a ij 


/ 


q Mi/2) 


3n, 


(ta r piq ) dq 


q i-0/2) 


ah <*bj - *m> ; ] * <ybj - ybi> m j 

1 4 




(Al) 


a n = * - Z a ik 

k*l J 

The last relation follows from the Gaussian condition (ref. 10); that is, 



In r 


Pq 


dq = 


7T 


To evaluate the b^. coefficients Simpson's rule is used for the case i + j and 
closed form integration is used for the case i = j since the integrand is singular for 
i=j. The result is 


h i 

b. =-!ln 
1} 6 


r. 


i, j-(l/2) + 4 1x1 r ij + 111 r i,j+(l/2) 


i * j 


(A2) 


/ h y 

b.. = h . I In — - 1 

11 1 V 2 ) 


11 


R. 


1 = ' 7 £ F kZ ^ r ik l AA k l 

q k,l 


(A3) 


where £ is the sum for all the plastic cells in the region and AA k j is the area of 


k, l 


the cell with coordinates (x k , y^). 

The coefficients A^ k , B i; j k> C.. k , and D.^ are given as follows using Simpson's 


rule: 


A *’* ft : { fc )2 : (y j - ^ t (x i - x bk> 2 - frj - W 

ljk 6| 4 

r ij,k-(l/2) 


r ijk 


(x i - ^+i> 2 - <y, - W 21 


r ij,k+(l/2) 


L 


B ■■ j (*j ■ ? k )(7i - -Ifc) t 4 (* t - x bk )( yi - y bk ) - jfr.jHyj - t 

ijk q ' 




(A4) 


4 

4 

4 


r 

ij,k-(l/2) ijk 

r ij,k+(l/2) 

J (A5) 

r h k 

J 

H 

M 

w 

i 

K* 

1 

x i - «k + i 


/a e\ 

6 

r 2 r>2 

+ 

2 


(A6) 


L ij>k-(l/2) r ijk 

r ij,k+(l/2)_ 



n _ X i ' X k 

u m 

o 


(A7) 


(Xi-x k ) + (yj - 





19 ! 



/ 


REFERENCES 

1. Mendelson, Alexander: Plasticity; Theory and Application. The Macmillan Co. , 

1968. 

2. Mendelson, Alexander: Elastic- Plastic Torsion Problems for Strain- Hardening 

Materials. NASA TN D-4391, 1968. 

3. Herakovich, C. T. : Strain Hardening Torsion of Solid Bars as a Minimization 

Problem. Int. J, Mech. Sci., vol. 12, no. 11, Nov. 1970, pp. 985-995. 

4. Herakovich, C. T. : Plastic Strains Associated with Torsion of Noncircular Bars. 

J. Appl. Mech., vol. 39, no. 1, Mar. 1972, pp. 285-287. 

5. Mendelson, Alexander: Boundary-Integral Methods in Elasticity and Plasticity. 

NASA TN D-7418, 1972. 

6. Rzasnicki, Walter; Mendelson, Alexander; and Albers, Lynn U. : Application of 

Boundary Integral Method to Elastoplastic Analysis of V-Notched Beams. NASA 
TN D-7424, 1973. 

7. Sokolnikoff, Ivan S. : Mathematical Theory of Elasticity. 2nd ed. , McGraw-Hill 

Book Co. , Inc., 1956. 

8. Prager, W. : An Introduction to the Mathematical Theory of Plasticity. J. Appl. 

Phys. , vol. 18, no. 4, April, 1947, pp. 375-383. 

9. Huth, J. H. : A Note on Plastic Torsion. J. Appl, Mech., vol. 22, no. 3, Sept. 

1955, pp. 432-434. 

10. Jaswon, M. A. ; and Ponter, A. R. : An Integral Equation Solution of the Torsion 
Problem. Proc. Roy. Soc. (London), Ser. A, vol. 273, no. 1353, May 1963, 
pp. 237-246. 


13 



TABLE I. - COMPARISON OF VALUES OF DIMENSIONLESS WARPING FUNCTION ON BOUNDARY 
OF ELASTIC SQUARE PLATE WITH EXACT ELASTIC SOLUTION 


Boundary 

value, 

y 

Exact 

warping 

function, 

w 

Value of warping function by 
boundary integral method 

Boundary 
value , 

y 

Exact 

warping 

function, 

W 

Value of warping function by 
boundary integral method 

Intervals, n 

Intervals, n 

4 


12 

16 

4 

8 

12 

16 

0.03125 

0.01095 




0.01095 

0. 4313 

0.1424 




0. 1424 

. 04167 

.01459 



0. 01459 


.5417 

.1433 



0. 1433 


.0625 

.02185 


0.02184 



.5625 

.1448 


0. 1446 



.09375 

.03264 




. 03264 

.5938 

. 1461 




. 1460 

.1250 

. 04328 

0.04311 


. 04328 


. 6250 

.1461 

0. 1446 


. 1461 


.1563 

.05374 




. 05372 

. 6563 

.1449 




. 1449 

.1875 

. 06390 


.06386 



. 6875 

.1422 


. 1419 



.2083 

. 07050 



.07051 


. 7083 

. 1396 



. 1395 


.2188 

.07380 




. 07377 

.7188 

.1380 




. 1380 

.2813 

, 09235 




. 09234 

.7813 

.1244 




. 1244 

.2917 1 

. 09527 



.09525 


. 7917 

.1214 



, 1212 i 


.3125 j 

.1009 


.1009 



.8125 

.1148 


. 1140 



.3438 j 

.1090 




.1090 

.8438 

. 1029 




. 1027 

.3750 

. 1 1 65 

, H59 


.116-1 


.5750 

.08664 

. 08826 


.08811 


.4063 

.1233 




. 1232 

.9063 

.07166 




.07129 

.4375 1 

.1293 


.1292 



.9375 

.05169 


.05142 



. 4583 1 

. 1329 



.1329 


.9583 

. 03644 



.03621 


. 4688 i 

. 1346 




.1346 

.9688 

. 02808 




. 02796 


TABLE II. - COMPARISON OF ELASTIC SOLUTIONS FOR MAXIMUM 
DIMENSIONLESS SHEAR STRESS AND DIMENSIONLESS 
MOMENT FOR SQUARE BAR 



Exact solution 

Finite difference 
method 
(55 eqs. ) 

Boundary integral method 

Intervals, n 

4 

8 

Dimensionless moment, M* 
Dimensionless maximum 
shear stress, r 

Hi it X 

1. 125 
.6754 

1. 122 
.6725 

1. 128 
.6724 

1. 127 
.6747 
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TABLE III. - COMPARISON OP ELASTIC SOLUTIONS 


FOR DIMENSIONLESS x DIRECTIONAL SHEAR 

STRESS DISTRIBUTION r x 

[First number, exact; second, boundary integral 
method; third, finite difference method.] 


y 

X 


0 

0.2 

0.4 

0.6 

0.8 

1, 0 


Elastic solutions 

1. 0 

0. 675 

0.658 

— 

0.605 

0.507 

0.342 

0 


.675 

.658 

. 605 

.506 

.339 

1 


.671 

.654 

. 600 

.500 

. 330 

1 

0.8 

0.492 

0.476 

0. 427 

0. 338 

0. 198 

0 


.492 

.476 

.427 

.338 

. 199 

1 


.492 

. 476 

. 428 

.339 

.200 

T 

0.6 

0,339 

0.326 

0.287 

0.219 

0. 121 

0 


.339 

.326 

.287 

.220 

. 122 

1 


.340 

.327 

.288 

.220 

. 123 

* 

0.4 

0.212 

0.203 

0. 176 

0. 132 

0. 0714 

0 


, .212 

.203 

. 177 

. 132 

. 0717 

I 


.212 

.203 

, 177 

. 133 

.0720 

1 

0.2 

0. 101 

0. 0971 

0.0839 

0. 0623 

0.0333 : 

0 


. 101 

.0971 

.0840 

.0624 

. 0335 

1 


. 102 

.0973 

.0842 

.0626 

.0336 





TABLE IV- - SUMMARY OF RESULTS FOR TORSION 


OF SQUARE PRISMATIC BAR 


Linear 
strain- 
hardening 
parameter, 
m 1 

Dimensionless 
angle of twist 
per unit 
length, 

0 

Dimensionless 

moment, 

M* 

Maximum 
dimensionless 
shear stress, 

r niax 

Maximum 

dimensionless 

strain, 

yP 

p _ 'max 

nnax 

G o 

0 

2 

1.81 

0.751 

0.792 


3 

1.96 



1.70 


4 

2.02 



2.57 


5 

2.06 



3.42 


6 

2.08 



4.25 

0.05 

2 

1.83 

0. 784 
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Figure 2. - Region R, don ndary curve c, and 
geometric quantities entering into boundary 
integrals. 
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figure 6, - Plastic zone boundaries in quadrant of square cross section as func- 
tion of dimensionless angle of twist per unit length for strain-hardening pa- 
rameter, 0.1. 
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figure 7. - Boundary and interior rotation for computing co- 
efficients given in appendix. 
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